model{
 for(i in 1:n){# n is the number of obs

  xb1[i] <- alpha1 + beta1[1]*x1[i,1] + beta1[2]*x1[i,2] + beta1[3]*x1[i,3] 
  + beta1[4]*x1[i,4] + beta1[5]*x1[i,5] + beta1[6]*x1[i,6] + beta1[7]*x1[i,7] 
  + beta1[8]*x1[i,8] + beta1[9]*x1[i,9] + beta1[10]*x1[i,10] 
  + beta1[11]*x1[i,11] 
  y1[i] ~ dnorm(xb1[i], tau1)

  xb2[i] <- alpha2 + beta2[1]*x2[i,1] + beta2[2]*x2[i,2] + beta2[3]*x2[i,3] 
  + beta2[4]*x2[i,4] + beta2[5]*x2[i,5] + beta2[6]*x2[i,6] + beta2[7]*x2[i,7]  
  + beta2[8]*x2[i,8] + beta2[9]*x2[i,9] + beta2[10]*x2[i,10] 
  + beta2[11]*x2[i,11] + beta2[12]*x2[i,12]
  y2[i] ~ dnorm(xb2[i], tau2) 

 }

  for(j in 1:K1){
   beta1[j] ~ dnorm(0, 0.1)
  }
  for(j in 1:K2){
   beta2[j] ~ dnorm(0, 0.1)
  }

  alpha1 ~ dnorm(0, 0.1)
  alpha2 ~ dnorm(0, 0.1)

# rho is the covariance between equations
 rho ~ dunif(-1,1) 
 
# partially observed covariance matrix 
 b0[1] <- 0
 b0[2] <- 0 
 B0[1,1] <- 1
 B0[2,2] <- 1
 B0[1,2] <- rho
 B0[2,1] <- rho

 prec[1:2,1:2] <- inverse(B0[,])
 SIGMA[1:2] ~ dmnorm(b0, prec)

 sigma1 <- SIGMA[1]
 sigma2 <- SIGMA[2]
 tau1 <- exp(sigma1)
 tau2 <- exp(sigma2)
}